Preliminares
En este documento ajustaremos algunos modelos de regresión lineal a los datos sobre venta de vinos. Para ello, utilizamos el conjunto de datos que generamos tras la depuración, asegurando un conjunto de datos “limpios” y exentos de ciertos peligros.
En primer lugar leemos el archivo generado en la depuración. Nótese que aquí se podría “repetir” el proceso con distintas posibilidades de tratamiento previo de las variables. Tal vez probar con winsorize e imputaciones aleatorias y luego comparar con otro esquema como el de eliminición de outliers y missing o con conversión de outliers a missings e imputación por modelos multivariantes…Todo un mundo de posibilidades!!
import pandas as pd
import numpy as np
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from sklearn.linear_model import LinearRegression
# Leer datos depurados datosvinoDep
vinosDep = pd.read_csv('C:\\Users\\Guille\\Documents\\MineriaDatos_2022_23\\PARTE I_Depuracion y Regresiones\\Dia1_MDDepuracion\\DatosVinoDep_winsRand.csv', index_col=0)
# Descriptivo de comprobación
vinosDep.head()## ID Acidez AcidoCitrico Azucar ... Region prop_missings Beneficio Compra
## 0 2 0.16 -0.81 26.10 ... 1.0 7.692308 515 1
## 1 4 2.64 -0.88 14.80 ... 3.0 0.000000 585 1
## 2 8 0.29 -0.40 21.50 ... 1.0 0.000000 0 0
## 3 11 -1.22 0.34 1.40 ... 2.0 7.692308 775 1
## 4 12 0.27 1.05 11.25 ... 2.0 0.000000 596 1
##
## [5 rows x 17 columns]Solicitamos información sobre el tipo de variables del archivo a modo de comprobación y observamos que se han cambiado los tipos! Esto se produce en el proceso de escritura y lectura del csv… Habrá que cambiarlo.
vinosDep.info()## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 6365 entries, 0 to 6364
## Data columns (total 17 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 ID 6365 non-null int64
## 1 Acidez 6365 non-null float64
## 2 AcidoCitrico 6365 non-null float64
## 3 Azucar 6365 non-null float64
## 4 CloruroSodico 6365 non-null float64
## 5 Densidad 6365 non-null float64
## 6 pH 6365 non-null float64
## 7 Sulfatos 6365 non-null float64
## 8 Alcohol 6365 non-null float64
## 9 CalifProductor 6365 non-null int64
## 10 PrecioBotella 6365 non-null float64
## 11 Etiqueta 6365 non-null object
## 12 Clasificacion 6365 non-null object
## 13 Region 6365 non-null float64
## 14 prop_missings 6365 non-null float64
## 15 Beneficio 6365 non-null int64
## 16 Compra 6365 non-null int64
## dtypes: float64(11), int64(4), object(2)
## memory usage: 895.1+ KB# Descripción de los datos
vinosDep.describe(include='all')## ID Acidez ... Beneficio Compra
## count 6365.000000 6365.000000 ... 6365.000000 6365.000000
## unique NaN NaN ... NaN NaN
## top NaN NaN ... NaN NaN
## freq NaN NaN ... NaN NaN
## mean 8010.702278 0.330902 ... 452.380204 0.785232
## std 4654.939139 0.768507 ... 308.380542 0.410694
## min 2.000000 -2.050000 ... 0.000000 0.000000
## 25% 3980.000000 0.130000 ... 236.000000 1.000000
## 50% 8065.000000 0.280000 ... 480.000000 1.000000
## 75% 12027.000000 0.650000 ... 671.000000 1.000000
## max 16128.000000 2.680000 ... 1568.000000 1.000000
##
## [11 rows x 17 columns]# Lista de columnas con menos de 10 valores distintos. Potenciales factores!
to_factor = list(vinosDep.loc[:,vinosDep.nunique() < 10]);
# Podemos cambiar el tipo de todas ellas a factor de una vez
vinosDep[to_factor] = vinosDep[to_factor].astype('category')
# Ordenmaos categorías de los fatcores de interés
vinosDep["Etiqueta"] = vinosDep["Etiqueta"].cat.reorder_categories(['MM','M','R','B','MB'])
vinosDep["Clasificacion"] = vinosDep["Clasificacion"].cat.reorder_categories(['Desc','*','**','***','****'])Estudio descriptivo de relaciones con la respuesta
En este apartado intentaremos descubrir a priori las relaciones marginales de las variables con la variable objetivo para hacernos una idea de cuales de ellas serán potencialmente influyentes en los modelos de regresión que ajustemos.
import seaborn as sns
# g = sns.PairGrid(vinosDep.iloc[:,6:])
# g.map_diag(sns.histplot)
# g.map_lower(sns.kdeplot)
# g.map_upper(sns.scatterplot)
# g.add_legend()
#
# del gPodemos generar este tipo de gráficos de rejilla, donde se pueden visualizar las relaciones a pares de todas las variables del archivo. En este caso, dado el caracter de los gráficos, solo será de aplicación a variables de tipo numérico.
No tiene muy buena pinta ya que las variables no parecen generar un patrón especial frente a la varuiable objetivo…se intuye poco valor predictivo de las variables para la regresión lineal.
Por otro lado, llama la antención la gran cantidad de 0 que tiene la variable objetivo continua Beneficio…Distribución muy compleja de modelizar..Tal vez habría que replantear el análisis y generar un modelo de regresión lineal para estimar el beneficio de los vinos que presentan algún beneficio.. Teniendo así una distribución de la variable respuesta más llevadera.
En este punto, podríamos pensar en una doble vía de actuación para el modelo predictivo final que generemos. Por una parte, un modelo de clasificación previa que estime la probabilidad de ser un vino con Compra = 1 (todo vino que tiene beneficio > 0) y posterioremente, para los clasificados como 1 en compra, un modelo de regresión lineal que estime el beneficio esperado. De alguna forma, estamos planteando un ensamble de modelos. Solo por plantear alternativas.
Variables de control
No es mala idea generar un par de variables de “control” para la evaluación de los efectos de los predictores frente a la respuesta. La idea es la siguiente: si generamos variables en el más estricto sentido aleatorio (por ejemplo siguiendo una distribución uniforme[0,1]) cualquier relación que estas presenten con la variable respuesta serán debidas puramente al azar, con lo que se pueden considerar relaciones espurias, es decir, falsas.
Por tanto, ya sea en la inspección preliminar de relaciones con la respuesta mediante correlación (relación lineal, válido para continua-continua) o VCramer (asociación en tablas de contingencia, válido para cruce de variables categóricas/nominales o continuas tramificadas) o bien en los propios modelos de regresión, las variables que presenten una menor relación con la respuesta que las variables de control, tendrán una sombra de sospecha sobre la veracidad de esa relación y probablemente serán descartadas, al menos en su estado original (siempre se pueden tratar de transformar, tramificar etc)
vinosDep['aleatorio'] = np.random.uniform(0,1,size=vinosDep.shape[0])
vinosDep['aleatorio2'] = np.random.uniform(0,1,size=vinosDep.shape[0])
vinosDep.head()## ID Acidez AcidoCitrico Azucar ... Beneficio Compra aleatorio aleatorio2
## 0 2 0.16 -0.81 26.10 ... 515 1 0.097808 0.593139
## 1 4 2.64 -0.88 14.80 ... 585 1 0.366505 0.062341
## 2 8 0.29 -0.40 21.50 ... 0 0 0.502511 0.794030
## 3 11 -1.22 0.34 1.40 ... 775 1 0.928091 0.538947
## 4 12 0.27 1.05 11.25 ... 596 1 0.684028 0.511079
##
## [5 rows x 19 columns]Interesante función para obetener un reporte descriptivo del archivo completo que evalúa a nivel univariante y bivariante las variables presentes en el archivo. También da información sobre valores extremos y missings, por lo que la podríamos utilizar en nuestra etapa de depuración de los datos.
# pip install pandas_profiling
import pandas_profiling
#pandas_profiling.ProfileReport(vinosDep)Muchísima información disponible en este reporte automático. A nivel univariante podemos comprobar distribuciones de variables continuas y representatividad de las categorías de los factores. A nivel bivariado, se evalúan las relaciones a pares de distintas formas. Especialmente interesante el gráfico de correlación Auto por su filosofía. Genera la correlación de spearman para continua-continua y v de cramer para categórica-categórica y para continua-categórica previamente tramificada la primera. Va muy en nuestra línea de pensamiento.
Con ese 21% de valores nulos, vamos a decidir partir el archivo para considerar solamente los vinos de Compra = 1 y ajustaremos el modelo de regresión sobre esta submuestra. Es importante que estas decisiones sean consensuadas y se especifiquen claramente en el informe de análisis, trasladándose a la fase de interpretación de resultados.
Creamos el archivo vinosCompra a este efecto.
vinosCompra = vinosDep.loc[vinosDep.Compra == 1].drop(['Compra'],axis=1)
vinosCompra.info()## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 4998 entries, 0 to 6364
## Data columns (total 18 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 ID 4998 non-null int64
## 1 Acidez 4998 non-null float64
## 2 AcidoCitrico 4998 non-null float64
## 3 Azucar 4998 non-null float64
## 4 CloruroSodico 4998 non-null float64
## 5 Densidad 4998 non-null float64
## 6 pH 4998 non-null float64
## 7 Sulfatos 4998 non-null float64
## 8 Alcohol 4998 non-null float64
## 9 CalifProductor 4998 non-null int64
## 10 PrecioBotella 4998 non-null float64
## 11 Etiqueta 4998 non-null category
## 12 Clasificacion 4998 non-null category
## 13 Region 4998 non-null category
## 14 prop_missings 4998 non-null category
## 15 Beneficio 4998 non-null int64
## 16 aleatorio 4998 non-null float64
## 17 aleatorio2 4998 non-null float64
## dtypes: category(4), float64(11), int64(3)
## memory usage: 606.0 KB#pandas_profiling.ProfileReport(vinosCompra)En el reporte para el nuevo archivo, volvemos a tener alta correlación (en un sentido amplio) de las variables Etiqueta y Clasificación con la variable objetivo, apareciendo algo coloreada la variable Alcohol. Pocas relaciones destacables a parte de esto.
Ya que vamos a hacer cosas como evaluación de las relaciones entre los predictores y la respuesta o creación masiva de transformaciones para conseguir linealidad, lo mejor es separar las respuestas y quedarnos con el input depurado, de esta forma podemos aplicar una misma función a todo el conjunto sin peligro de transformar las respuestas y cosas raras que puedan suceder.
# Eliminar variable objetivo continua
varObjCont = vinosCompra.Beneficio
imputCompra = vinosCompra.drop(['ID','Beneficio'],axis=1)
varObjCont.info()## <class 'pandas.core.series.Series'>
## Int64Index: 4998 entries, 0 to 6364
## Series name: Beneficio
## Non-Null Count Dtype
## -------------- -----
## 4998 non-null int64
## dtypes: int64(1)
## memory usage: 78.1 KBNos encantaría tener un procedimiento que evaluara las relaciones de todas las variables del input frente a la variable objetivo (ya sea continua o binaria) y que de alguna forma ordenara esos efectos por intensiadad de asociación para generar un ranking tentativo a priori de las vaiables más interesantes a nivel individual (es importante destacar que estas relaciones solamente cuentan el efecto marginal de una variable con la objetivo sin considerar las posibles interacciones existentes entre variables que pueden modificar en modelo estas asociaciones).
Vamos a definir algunas funciones que nos faciliten el proceso de búsqueda, automatizando lo más posible para que luego sea aplicar y evaluar!
Primero nos generamos la función cramers_v, cuyo objetivo es calcular el valor de v de cramer para la asociación entre dos variables cualesquiera (se entiende que una de ellas será una variable objetivo…pero podría funcionar entre predictores de igual forma..) Molaría que la función distinguiera el caso de variable continua y categórica puesto que la asociación mediante vCramer se entiende en términos de cruce de dos factores y, por tanto, si la variable es continua, debe existir un proceso previo de tramificación o discretización de la variable. Entonces, programamos el proceso de tal manera que si encuentra columnas numéricas las tramifique (en 5 tramos por ejemplo) y sino simplemente calcule el valor cramer que se extrae de la tabla de contingencia de los dos factores introducidos.
La filosofía de programación es la misma que en otras funciones definidas anteriormente. Como realmente querremos aplicarla sobre el input completo, estamos pensando en utilizar un apply para que aplique un mismo proceso a todas las columnas que vaya encontrando, manteniendo la variable objetivo. Por ello, es importante separar esta objetivo del conjunto de predictores así como lo hacíamos en imputaciones y transformación de missings.
# Librería estadística!
import scipy.stats as stats
# Función para calcular VCramer (dos nominales de entrada!)
def cramers_v(var1, varObj):
if not var1.dtypes == 'category':
bins = min(5,var1.value_counts().count())
var1 = pd.cut(var1, bins = 5)
if not varObj.dtypes == 'category': #np.issubdtype(varObj, np.number):
bins = min(5,varObj.value_counts().count())
varObj = pd.cut(varObj, bins = 5)
data = pd.crosstab(var1, varObj).values
vCramer = stats.contingency.association(data, method = 'cramer')
return vCramer
# Ejemplo uso univariante
#cramers_v(vinosCompra['Etiqueta'],vinosCompra['Beneficio'])
# Aplicar la función al input completo contra la objetivo
tablaCramer = pd.DataFrame(imputCompra.apply(lambda x: cramers_v(x,varObjCont)),columns=['VCramer'])
# Obtener el gráfico de importancia de las variables frente a la objetivo continua según vcramer
import plotly.express as px
px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones frente al Beneficio').update_yaxes(categoryorder="total ascending")Ya tenemos nuestro ranking de variables según cramer para todos los predictores en relación a la variable objetivo continua, Beneficio. Como ya intuíamos, las variables categóricas Etiqueta y Clasificación son las que mayor relación presentan, es decir, las que mayor patrón de asociación pueden aportar para el futuro modelo.
Entendiendo el test Chi-cuadrado
Vamos a realizar un análisis de asociación de una tabla de contingencia basada en Chi-cuadrado. Para ilustrar esto, nos permitimos el lujo de adelantar algo que observaremos posteriormente y que podría afectar a nuestro modelo por la presencia de interacción entre las variables más relevantes Clasificación y Etiqueta.
# Opción para calcular las medias de Beneficio en los cruces de Clasificación y etiqueta
pd.pivot_table(vinosCompra,index =['Etiqueta'],columns=['Clasificacion'],values=['Beneficio'],aggfunc=[np.mean])
# Tabla de contingencia con la cuenta de frecuencia de observaciones en cada cruce## mean
## Beneficio
## Clasificacion Desc * ** *** ****
## Etiqueta
## MM 227.714286 292.464646 349.676471 352.909091 NaN
## M 365.428571 403.221445 449.474178 523.120968 635.846154
## R 507.620553 536.413989 581.393564 636.366733 725.234694
## B 706.515625 635.623457 701.191919 775.459016 857.936416
## MB 900.750000 830.000000 826.568182 857.851351 929.973684pd.crosstab(vinosCompra.Etiqueta,vinosCompra.Clasificacion)## Clasificacion Desc * ** *** ****
## Etiqueta
## MM 63 99 34 11 0
## M 273 429 426 124 13
## R 253 529 808 499 98
## B 64 162 396 366 173
## MB 12 10 44 74 38Vayamos un poco más allá por un momento y nos planteamos si habrá un efecto de interacción entre estas dos variable relevantes. Que demonios queremos decir con esto? Si hay interacción, es posible que el efecto conjunto de las dos variable sobre la respuesta no sea uniforme o bien no sea el resultado de la adición simple de ambos efectos. De esta forma, nos planteamos si es igualmente probable encontrar un vino con Clasificación máxima y la peor Etiqueta, que un vino con Clasificación máxima y la mejor de las Etiquetas… La lógica puede inducir a pensar que debería existir cierto patrón de asociación..
# Probamos el test chi cuadrado para ver las salidas
stats.chi2_contingency(pd.crosstab(vinosCompra.Etiqueta,vinosCompra.Clasificacion).values)
#pd.crosstab(data_trans.Beneficio_cut,vinosCompra.Clasificacion,margins=True,normalize='columns')## (842.5114902984811, 5.340776470642584e-169, 16, array([[ 27.54201681, 50.90096038, 70.7394958 , 44.48139256,
## 13.33613445],
## [168.31232493, 311.06142457, 432.29691877, 271.83073229,
## 81.49859944],
## [290.98739496, 537.77971188, 747.37815126, 469.95558223,
## 140.89915966],
## [154.47478992, 285.4879952 , 396.75630252, 249.48259304,
## 74.79831933],
## [ 23.68347339, 43.76990796, 60.82913165, 38.24969988,
## 11.46778711]]))- Primer valor: Valor del estadístico Chi-cuadrado
- Segundo valor: p-valor del contraste de hipótesis (recordemos H0: no existe asociación)
- Tercer valor: Grados de libertad de la distribución de contraste (filas-1)*(columnas-1)
- Cuarto valor: Matriz de valores esperados bajo H0 (según las distribuciones marginales de las variables, cuan probable es encontrar casilla 1,1… p(x=1)*p(y=1)
El valor de chi es muy alto, por lo que el p-valor es muy pequeño, prácticamente 0. Por lo que se evidencia una asociación significativa entre estas dos variables, o bien, la matriz de valores observados es significativamente distinta a la matriz de valores esperados bajo independencia.
Precaución!! Hay una cosa importante en el test Chi-cuadrado. Los resultados no son robustos si alguna de las casillas tiene frecuencia esperada menor que 5. Si se da este caso, sería necesario recodificar alguna de las variables del cruce para obtener tamaño muestral suficiente en las marginales. En nuestro caso, esto no llega a suceder (ver matriz de frecuencias esperadas), sin embargo tenemos una casilla con frecuencia observada 0, lo cual puede estar tirando mucho del valor de la chi y de alguna forma distorsionando sus resultados. Por ello, no es mala práctica unir alguna categoría para evitar esto.
Veamos los cruces de estas categóricas con Beneficio.
#sns.boxplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
print(sns.violinplot(x='Etiqueta',y='Beneficio',data=vinosCompra,palette='viridis'))
#sns.stripplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
#sns.swarmplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')## AxesSubplot(0.125,0.11;0.775x0.77)#sns.boxplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
#sns.violinplot(x='Etiqueta',y='Beneficio',data=vinosCompra,palette='viridis')
sns.stripplot(x='Clasificacion',y='Beneficio',data=vinosDep,palette='viridis')
#sns.swarmplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
Evaluación de la asociación categorica-numérica mediante ANOVA
Los contrastes estadísticos pueden resultar de utilidad para contatar la significación de las relaciones visualizadas a nivel gráfico. De esta forma, si tenemos una factor y queremos cuantificar esta significación con respecto a la variable objetivo continua, no es descabellado pensar en un análisis de la varianza con un factor para valorar si las distribuciones en los subgrupos difieren de la distribución global.
El Anova realiza la descomposición de la suma de cuadrados en las componentes explicada y residual, de tal forma que imputamos una parte de la variabilidad de la variable continua al cambio en los niveles del factor y la parte que no se puede explicar solamente por ese cambio, queda como variabilidad no explicada. La idea es que la variabilidad explicada sea mucho mayor que la residual.
Para la variable Región: Se observa que el p-valor es superior a 0,05 por lo que la asociación no es significativa al 95% de nivel de confianza. Esta variable no influye mucho en el Beneficio.
# get ANOVA table as R like output
import statsmodels.api as sm
from statsmodels.formula.api import ols
# ANOVA región
model = ols('Beneficio ~ C(Region)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table## sum_sq df F PR(>F)
## C(Region) 2.362961e+05 2.0 2.372917 0.093313
## Residual 2.487021e+08 4995.0 NaN NaNPara la variable Etiqueta: Se observa que el p-valor es 0 por lo que la asociación es significativa a cualquier nivel de confianza. Esta variable influye mucho en el Beneficio.
# ANOVA Etiqueta
model = ols('Beneficio ~ C(Etiqueta)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table## sum_sq df F PR(>F)
## C(Etiqueta) 9.240747e+07 4.0 736.899964 0.0
## Residual 1.565309e+08 4993.0 NaN NaNANOVA de doble vía. Puede resultar útil para la evaluación de efectos conjuntos de factores sobre la objetivo continua. En este caso, evaluamos la interacción de la que ya sospechábamos entre Clasificación y Etiqueta.
# ANOVA con dos factores e interacción
model = ols('Beneficio ~ C(Etiqueta) + C(Clasificacion) + C(Etiqueta):C(Clasificacion)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table## sum_sq df F PR(>F)
## C(Etiqueta) 7.500830e+07 4.0 659.687718 0.000000e+00
## C(Clasificacion) 4.971909e+06 4.0 43.727257 4.012464e-36
## C(Etiqueta):C(Clasificacion) 1.016574e+07 16.0 22.351569 7.635523e-64
## Residual 1.413894e+08 4974.0 NaN NaNSe constata la significación estadística de todos los efectos presentes en el modelo. Podríamos decir que existe influencia marginal de ambos factores sobre el Beneficio pero también una influencia conjunta de ambos! Así, el beneficio esperado puede depender de la combinación de Clasificación y Etiqueta del vino.
Para ilustrar este efecto, podemos recurrir al gráfico de interacciones entre dos factores y la variable continua. Aquí, detalle de programación y es que no acepta el tipo category (con su orden) por lo que habrá que convertir las columnas en tipo ‘object’ (pero solo para esto!)
from statsmodels.graphics.factorplots import interaction_plot
import matplotlib.pyplot as plt
# Gráfico de interacción de las componentes del ANOVA
fig = interaction_plot(x=vinosCompra.Etiqueta.astype('object'), trace=vinosCompra['Clasificacion'].astype('object'), response=vinosCompra['Beneficio'])
plt.show()
Se evidencia interacción entre los factores cuando las líneas se cortan en algún punto. Esto quier decir que la evolución del Beneficio a lo largo de las Etiquetas no se produce de la misma forma en todas los niveles de Clasificación.
Se puede observar esa falta de dato en la combinación MM - ****
Vamos a hacer la prueba de recategorizar alguno de los niveles implicados. Para decidir cual de ellos, recurriremos a 1) frecuencia de la categoría (reporte) y 2) Similitud con la categoría adyacente (boxplot-violin). De esta forma decidimos que parece más conveniente unir Etiqueta en sus categorías M y MM.
# Unimos la categoría minoritaria MM con M en Etiqueta
vinosCompra['Etiqueta_r'] = vinosCompra.Etiqueta.replace(['M','MM'],'M-MM',inplace=False)pd.crosstab(vinosCompra.Etiqueta_r,vinosCompra.Clasificacion)## Clasificacion Desc * ** *** ****
## Etiqueta_r
## M-MM 336 528 460 135 13
## R 253 529 808 499 98
## B 64 162 396 366 173
## MB 12 10 44 74 38# Probamos el test chi cuadrado para ver las salidas
stats.chi2_contingency(pd.crosstab(vinosCompra.Etiqueta_r,vinosCompra.Clasificacion).values)## (800.5770302943753, 1.244618385605778e-163, 12, array([[195.85434174, 361.96238495, 503.03641457, 316.31212485,
## 94.83473389],
## [290.98739496, 537.77971188, 747.37815126, 469.95558223,
## 140.89915966],
## [154.47478992, 285.4879952 , 396.75630252, 249.48259304,
## 74.79831933],
## [ 23.68347339, 43.76990796, 60.82913165, 38.24969988,
## 11.46778711]]))# Gráfico de interacción de las componentes del ANOVA
fig = interaction_plot(x=vinosCompra.Etiqueta_r.astype('object'), trace=vinosCompra['Clasificacion'].astype('object'), response=vinosCompra['Beneficio'])
plt.show()
Sigue evidenciándose interacción entre los factores, especialmente en la categoría desconocido.
Transformaciones de las variables continuas
El principal objetivo de las transformaciones de las variables continuas es conseguir linealidad frente a la variable objetivo. De esta forma se prueban las transformaciones típicas (log, exp, potencias y raíces) y se escoge aquella que mayor coeficiente de correlación o valor V de Cramer presenta con la respuesta continua (o binaria en su caso).
Definamos pues una función que nos facilite este proceso, de tal manera que aplique sobre una columna genérica y una variable objetivo e implemente la posibilidad de actuar en base al coeficiente de correlación (relación lineal de continua-continua) o al valor V de cramer (asociación general entre categórica-categórcia).
La función va a hacer una traslación a valores positivos de la variable que encuentre (predictor continuo siempre!! no tiene sentido el log de una categórica!) y seguidamente va a generar las posibles transfromaciones típicas (log, exp, potencias y raíces 2 y 4) para aplicar estos criterios sobre todas las posibilidades y escoger la que mejor valor obtenga. Vamos a ello!
Aquí hay un walkaround para conseguir que los nombres de las transfromadas se conserven en el dataset de salida que es fruto de un apply al input continuo. Seguro se puede mejorar!
# Busca la transformación de variables input de intervalo que maximiza la VCramer o la correlación tipo Pearson con la objetivo
def mejorTransf (vv,target, name=False, tipo = 'cramer', graf=False):
# Traslación a valores positivos de la variable (sino falla log y las raíces!)
vv = vv + abs(min(vv))+0.0001
# Definimos y calculamos las tranformacione típicas
transf = pd.DataFrame({vv.name + '_ident': vv, vv.name + '_log': np.log(vv), vv.name + '_exp': np.exp(vv),
vv.name + '_sqr': np.square(vv), vv.name + '_cuarta': vv**4, vv.name + '_raiz4': vv**(1/4)})
# Distinguimos caso cramer o caso correlación
if tipo == 'cramer':
# Aplicar la función cramers_v a cada trasnformación frente a la respuesta
tablaCramer = pd.DataFrame(transf.apply(lambda x: cramers_v(x,target)),columns=['VCramer'])
# Si queremos gráfico, muestra comparativa entre las posibilidades
if graf: px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones frente a ' + target.name).update_yaxes(categoryorder="total ascending").show()
# Identificar mejor transfromación
best = tablaCramer.query('VCramer == VCramer.max()').index
ser = transf[best].squeeze()
if tipo == 'cor':
# Aplicar coeficiente de correlacion a cada trasnformación frente a la respuesta
tablaCorr = pd.DataFrame(transf.apply(lambda x: np.corrcoef(x,target)[0,1]),columns=['Corr'])
# Si queremos gráfico, muestra comparativa entre las posibilidades
if graf: px.bar(tablaCorr,x=tablaCorr.Corr,title='Relaciones frente a ' + target.name).update_yaxes(categoryorder="total ascending").show()
# identificar mejor transfromación
best = tablaCorr.query('Corr.abs() == Corr.abs().max()').index
ser = transf[best].squeeze()
# Aquí distingue si se devuelve la variable transfromada o solamente el nombre de la transfromacion
if name:
return(ser.name)
else:
return(ser)
# Ejemplo de uso univariante
tr = mejorTransf(vinosCompra.Azucar,varObjCont, tipo='cor')
tr.head(),vinosCompra.Azucar.head()## (0 15055.31454
## 1 12409.98228
## 3 9604.01960
## 4 11631.64407
## 5 9682.57968
## Name: Azucar_sqr, dtype: float64, 0 26.10
## 1 14.80
## 3 1.40
## 4 11.25
## 5 1.80
## Name: Azucar, dtype: float64)#tr=pd.DataFrame()
#for i in range(len(imputCompra.select_dtypes(include=np.number))):
# tr.iloc[:,i] = mejorTransfVCramer(imputCompra[column],varObjCont)imputCompra.info()## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 4998 entries, 0 to 6364
## Data columns (total 16 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 Acidez 4998 non-null float64
## 1 AcidoCitrico 4998 non-null float64
## 2 Azucar 4998 non-null float64
## 3 CloruroSodico 4998 non-null float64
## 4 Densidad 4998 non-null float64
## 5 pH 4998 non-null float64
## 6 Sulfatos 4998 non-null float64
## 7 Alcohol 4998 non-null float64
## 8 CalifProductor 4998 non-null int64
## 9 PrecioBotella 4998 non-null float64
## 10 Etiqueta 4998 non-null category
## 11 Clasificacion 4998 non-null category
## 12 Region 4998 non-null category
## 13 prop_missings 4998 non-null category
## 14 aleatorio 4998 non-null float64
## 15 aleatorio2 4998 non-null float64
## dtypes: category(4), float64(11), int64(1)
## memory usage: 656.9 KB# Aplicar a las variables continuas la mejor transfromación según correlacion frente a varObjCont
transf_cor = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont, tipo='cor'))
# Pedir los nombres de las transromadas
transf_cor_names = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont,tipo='cor', name=True))
# Asignar nombres a las columnas de salida del proceso
transf_cor.columns = transf_cor_names.values
transf_cor## Acidez_sqr AcidoCitrico_exp ... aleatorio_log aleatorio2_cuarta
## 0 4.884542 4.179117 ... -2.323380 0.123920
## 1 21.997038 3.896583 ... -1.003377 0.000015
## 3 0.689066 13.198458 ... -0.074482 0.084480
## 4 5.382864 26.845548 ... -0.379560 0.068321
## 5 3.349266 13.875157 ... -2.186242 0.007613
## ... ... ... ... ... ...
## 6357 4.973346 13.198458 ... -0.119273 0.018796
## 6361 7.508148 10.278969 ... -2.864621 0.004033
## 6362 5.664876 3.669664 ... -1.735340 0.002771
## 6363 4.928844 3.127081 ... -0.338840 0.207337
## 6364 5.617374 10.592011 ... -3.952220 0.086490
##
## [4998 rows x 12 columns]transf_cor## Acidez_sqr AcidoCitrico_exp ... aleatorio_log aleatorio2_cuarta
## 0 4.884542 4.179117 ... -2.323380 0.123920
## 1 21.997038 3.896583 ... -1.003377 0.000015
## 3 0.689066 13.198458 ... -0.074482 0.084480
## 4 5.382864 26.845548 ... -0.379560 0.068321
## 5 3.349266 13.875157 ... -2.186242 0.007613
## ... ... ... ... ... ...
## 6357 4.973346 13.198458 ... -0.119273 0.018796
## 6361 7.508148 10.278969 ... -2.864621 0.004033
## 6362 5.664876 3.669664 ... -1.735340 0.002771
## 6363 4.928844 3.127081 ... -0.338840 0.207337
## 6364 5.617374 10.592011 ... -3.952220 0.086490
##
## [4998 rows x 12 columns]# Aplicar a las variables continuas la mejor transfromación según cramer frente a varObjCont
transf_cramer = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont, tipo='cramer'))
transf_cramer_names = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont,tipo='cramer', name=True))
transf_cramer.columns = transf_cramer_names.valuestransf_cramer## Acidez_sqr AcidoCitrico_cuarta ... aleatorio_cuarta aleatorio2_cuarta
## 0 4.884542 4.182786 ... 9.201856e-05 0.123920
## 1 21.997038 3.422026 ... 1.806986e-02 0.000015
## 3 0.689066 44.314531 ... 7.423561e-01 0.084480
## 4 5.382864 117.175386 ... 2.190970e-01 0.068321
## 5 3.349266 47.850783 ... 1.592608e-04 0.007613
## ... ... ... ... ... ...
## 6357 4.973346 44.314531 ... 6.205854e-01 0.018796
## 6361 7.508148 29.478015 ... 1.055951e-05 0.004033
## 6362 5.664876 2.856979 ... 9.669522e-04 0.002771
## 6363 4.928844 1.689553 ... 2.578546e-01 0.207337
## 6364 5.617374 31.025702 ... 1.362358e-07 0.086490
##
## [4998 rows x 12 columns]# Input con transfromadas según correlación
imput_transf_cor = imputCompra.join(transf_cor)
# Input con transfromadas según cramer
imput_transf_cramer = imputCompra.join(transf_cramer)imput_transf_cramer.info()## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 4998 entries, 0 to 6364
## Data columns (total 28 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 Acidez 4998 non-null float64
## 1 AcidoCitrico 4998 non-null float64
## 2 Azucar 4998 non-null float64
## 3 CloruroSodico 4998 non-null float64
## 4 Densidad 4998 non-null float64
## 5 pH 4998 non-null float64
## 6 Sulfatos 4998 non-null float64
## 7 Alcohol 4998 non-null float64
## 8 CalifProductor 4998 non-null int64
## 9 PrecioBotella 4998 non-null float64
## 10 Etiqueta 4998 non-null category
## 11 Clasificacion 4998 non-null category
## 12 Region 4998 non-null category
## 13 prop_missings 4998 non-null category
## 14 aleatorio 4998 non-null float64
## 15 aleatorio2 4998 non-null float64
## 16 Acidez_sqr 4998 non-null float64
## 17 AcidoCitrico_cuarta 4998 non-null float64
## 18 Azucar_raiz4 4998 non-null float64
## 19 CloruroSodico_cuarta 4998 non-null float64
## 20 Densidad_log 4998 non-null float64
## 21 pH_cuarta 4998 non-null float64
## 22 Sulfatos_sqr 4998 non-null float64
## 23 Alcohol_ident 4998 non-null float64
## 24 CalifProductor_log 4998 non-null float64
## 25 PrecioBotella_sqr 4998 non-null float64
## 26 aleatorio_cuarta 4998 non-null float64
## 27 aleatorio2_cuarta 4998 non-null float64
## dtypes: category(4), float64(23), int64(1)
## memory usage: 1.1 MBPodemos comprobar si las transformaciones han aumentado el valor de correlación lineal con la respuesta.
En el gráfico de correlaciones se intuye la baja relación lineal presente entre las variables del archivo 2 a 2. Esto nos hace pensar que las variables continuas, a diferencia de las categóricas, no van a tener demasiada influencia en los modelos de regresión frente a la respuesta (véase el poco color que presenta la fila 1). Por otro lado, podemos estar respirar con cierta tranquilidad ante el hecho de la baja relación entre los predictores (véase la ausencia de color en todo el gráfico), cosa que evitará problemas de colinealidad en los modelos.
Es evidente que, las relaciones de cada variable con su transformada han de ser muy elevadas, ya que hay un mecanismo que genera una en función exacta de la otra, por lo que esto es normal.
Nota: La principal precaución que hay que tener es no considerar un modelo completo con todas al mismo tiempo puesto que se pueden generar los problemas de colinealidad. Solamente utilizaremos el set completo de variables cuando hagamos un proceso de selección automática de variables, proceso en el cual se elegirán las que más R2 aporten al modelo.
# Matriz de correlaciones incluyendo transformaciones
corr = imput_transf_cor.join(varObjCont).corr()
corr.style.background_gradient(cmap='coolwarm').format(precision=3)
#sns.heatmap(corr, annot=True)| Acidez | AcidoCitrico | Azucar | CloruroSodico | Densidad | pH | Sulfatos | Alcohol | CalifProductor | PrecioBotella | aleatorio | aleatorio2 | Acidez_sqr | AcidoCitrico_exp | Azucar_sqr | CloruroSodico_log | Densidad_log | pH_log | Sulfatos_exp | Alcohol_raiz4 | CalifProductor_raiz4 | PrecioBotella_sqr | aleatorio_log | aleatorio2_cuarta | Beneficio | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Acidez | 1.000 | -0.023 | -0.022 | 0.006 | 0.007 | 0.015 | 0.027 | -0.003 | 0.009 | -0.002 | 0.005 | -0.008 | 0.961 | -0.014 | -0.030 | -0.007 | 0.007 | 0.017 | 0.010 | -0.007 | 0.003 | -0.004 | 0.013 | -0.004 | -0.044 |
| AcidoCitrico | -0.023 | 1.000 | 0.007 | -0.017 | -0.014 | -0.010 | 0.002 | 0.003 | 0.061 | -0.013 | 0.015 | -0.005 | -0.029 | 0.796 | 0.005 | -0.006 | -0.014 | -0.011 | 0.011 | 0.006 | 0.058 | -0.015 | 0.021 | -0.010 | 0.002 |
| Azucar | -0.022 | 0.007 | 1.000 | -0.010 | 0.013 | 0.003 | -0.001 | -0.028 | -0.003 | -0.020 | -0.011 | 0.014 | -0.023 | -0.001 | 0.956 | 0.004 | 0.013 | 0.002 | -0.009 | -0.032 | 0.005 | -0.018 | -0.016 | 0.004 | -0.013 |
| CloruroSodico | 0.006 | -0.017 | -0.010 | 1.000 | 0.030 | -0.016 | -0.004 | -0.019 | 0.014 | -0.007 | 0.004 | 0.004 | 0.005 | -0.016 | -0.014 | 0.628 | 0.030 | -0.018 | -0.009 | -0.018 | 0.023 | -0.007 | 0.009 | -0.001 | -0.003 |
| Densidad | 0.007 | -0.014 | 0.013 | 0.030 | 1.000 | -0.016 | -0.012 | 0.000 | 0.035 | 0.008 | 0.008 | 0.015 | 0.006 | -0.013 | 0.011 | -0.002 | 1.000 | -0.015 | -0.012 | 0.007 | 0.031 | 0.010 | 0.010 | 0.018 | -0.015 |
| pH | 0.015 | -0.010 | 0.003 | -0.016 | -0.016 | 1.000 | 0.004 | -0.017 | -0.078 | -0.010 | -0.003 | 0.013 | 0.015 | 0.012 | 0.002 | -0.001 | -0.016 | 0.990 | 0.000 | -0.015 | -0.080 | -0.009 | 0.000 | 0.025 | 0.018 |
| Sulfatos | 0.027 | 0.002 | -0.001 | -0.004 | -0.012 | 0.004 | 1.000 | -0.003 | 0.019 | 0.004 | 0.006 | -0.011 | 0.023 | -0.003 | -0.006 | 0.001 | -0.012 | 0.006 | 0.698 | -0.004 | 0.015 | 0.005 | 0.012 | -0.013 | -0.005 |
| Alcohol | -0.003 | 0.003 | -0.028 | -0.019 | 0.000 | -0.017 | -0.003 | 1.000 | -0.025 | 0.015 | 0.002 | 0.018 | 0.001 | 0.008 | -0.028 | -0.011 | 0.001 | -0.021 | -0.010 | 0.955 | -0.033 | 0.012 | 0.005 | 0.018 | 0.082 |
| CalifProductor | 0.009 | 0.061 | -0.003 | 0.014 | 0.035 | -0.078 | 0.019 | -0.025 | 1.000 | 0.019 | 0.023 | -0.016 | 0.005 | 0.048 | -0.005 | 0.002 | 0.035 | -0.076 | 0.007 | -0.018 | 0.890 | 0.022 | 0.024 | -0.026 | -0.052 |
| PrecioBotella | -0.002 | -0.013 | -0.020 | -0.007 | 0.008 | -0.010 | 0.004 | 0.015 | 0.019 | 1.000 | 0.007 | -0.012 | -0.006 | -0.018 | -0.016 | 0.015 | 0.008 | -0.012 | -0.002 | 0.010 | 0.022 | 0.978 | 0.015 | -0.015 | 0.011 |
| aleatorio | 0.005 | 0.015 | -0.011 | 0.004 | 0.008 | -0.003 | 0.006 | 0.002 | 0.023 | 0.007 | 1.000 | 0.016 | 0.004 | 0.017 | -0.008 | -0.013 | 0.008 | -0.003 | -0.007 | 0.009 | 0.013 | 0.002 | 0.872 | 0.018 | -0.001 |
| aleatorio2 | -0.008 | -0.005 | 0.014 | 0.004 | 0.015 | 0.013 | -0.011 | 0.018 | -0.016 | -0.012 | 0.016 | 1.000 | -0.007 | 0.002 | 0.014 | -0.019 | 0.016 | 0.013 | -0.015 | 0.012 | -0.015 | -0.014 | 0.021 | 0.866 | 0.017 |
| Acidez_sqr | 0.961 | -0.029 | -0.023 | 0.005 | 0.006 | 0.015 | 0.023 | 0.001 | 0.005 | -0.006 | 0.004 | -0.007 | 1.000 | -0.018 | -0.030 | -0.012 | 0.006 | 0.017 | 0.006 | -0.004 | -0.001 | -0.008 | 0.011 | -0.001 | -0.044 |
| AcidoCitrico_exp | -0.014 | 0.796 | -0.001 | -0.016 | -0.013 | 0.012 | -0.003 | 0.008 | 0.048 | -0.018 | 0.017 | 0.002 | -0.018 | 1.000 | -0.004 | -0.005 | -0.012 | 0.012 | -0.000 | 0.011 | 0.049 | -0.017 | 0.025 | -0.005 | 0.014 |
| Azucar_sqr | -0.030 | 0.005 | 0.956 | -0.014 | 0.011 | 0.002 | -0.006 | -0.028 | -0.005 | -0.016 | -0.008 | 0.014 | -0.030 | -0.004 | 1.000 | 0.005 | 0.011 | 0.001 | -0.013 | -0.033 | 0.005 | -0.012 | -0.012 | 0.002 | -0.013 |
| CloruroSodico_log | -0.007 | -0.006 | 0.004 | 0.628 | -0.002 | -0.001 | 0.001 | -0.011 | 0.002 | 0.015 | -0.013 | -0.019 | -0.012 | -0.005 | 0.005 | 1.000 | -0.002 | -0.003 | 0.004 | -0.014 | 0.014 | 0.016 | -0.003 | -0.017 | 0.009 |
| Densidad_log | 0.007 | -0.014 | 0.013 | 0.030 | 1.000 | -0.016 | -0.012 | 0.001 | 0.035 | 0.008 | 0.008 | 0.016 | 0.006 | -0.012 | 0.011 | -0.002 | 1.000 | -0.015 | -0.012 | 0.008 | 0.031 | 0.010 | 0.010 | 0.019 | -0.015 |
| pH_log | 0.017 | -0.011 | 0.002 | -0.018 | -0.015 | 0.990 | 0.006 | -0.021 | -0.076 | -0.012 | -0.003 | 0.013 | 0.017 | 0.012 | 0.001 | -0.003 | -0.015 | 1.000 | 0.001 | -0.017 | -0.079 | -0.011 | -0.001 | 0.026 | 0.019 |
| Sulfatos_exp | 0.010 | 0.011 | -0.009 | -0.009 | -0.012 | 0.000 | 0.698 | -0.010 | 0.007 | -0.002 | -0.007 | -0.015 | 0.006 | -0.000 | -0.013 | 0.004 | -0.012 | 0.001 | 1.000 | -0.010 | 0.007 | 0.000 | -0.002 | -0.008 | -0.009 |
| Alcohol_raiz4 | -0.007 | 0.006 | -0.032 | -0.018 | 0.007 | -0.015 | -0.004 | 0.955 | -0.018 | 0.010 | 0.009 | 0.012 | -0.004 | 0.011 | -0.033 | -0.014 | 0.008 | -0.017 | -0.010 | 1.000 | -0.024 | 0.007 | 0.015 | 0.015 | 0.083 |
| CalifProductor_raiz4 | 0.003 | 0.058 | 0.005 | 0.023 | 0.031 | -0.080 | 0.015 | -0.033 | 0.890 | 0.022 | 0.013 | -0.015 | -0.001 | 0.049 | 0.005 | 0.014 | 0.031 | -0.079 | 0.007 | -0.024 | 1.000 | 0.023 | 0.015 | -0.021 | -0.053 |
| PrecioBotella_sqr | -0.004 | -0.015 | -0.018 | -0.007 | 0.010 | -0.009 | 0.005 | 0.012 | 0.022 | 0.978 | 0.002 | -0.014 | -0.008 | -0.017 | -0.012 | 0.016 | 0.010 | -0.011 | 0.000 | 0.007 | 0.023 | 1.000 | 0.010 | -0.019 | 0.012 |
| aleatorio_log | 0.013 | 0.021 | -0.016 | 0.009 | 0.010 | 0.000 | 0.012 | 0.005 | 0.024 | 0.015 | 0.872 | 0.021 | 0.011 | 0.025 | -0.012 | -0.003 | 0.010 | -0.001 | -0.002 | 0.015 | 0.015 | 0.010 | 1.000 | 0.027 | -0.004 |
| aleatorio2_cuarta | -0.004 | -0.010 | 0.004 | -0.001 | 0.018 | 0.025 | -0.013 | 0.018 | -0.026 | -0.015 | 0.018 | 0.866 | -0.001 | -0.005 | 0.002 | -0.017 | 0.019 | 0.026 | -0.008 | 0.015 | -0.021 | -0.019 | 0.027 | 1.000 | 0.027 |
| Beneficio | -0.044 | 0.002 | -0.013 | -0.003 | -0.015 | 0.018 | -0.005 | 0.082 | -0.052 | 0.011 | -0.001 | 0.017 | -0.044 | 0.014 | -0.013 | 0.009 | -0.015 | 0.019 | -0.009 | 0.083 | -0.053 | 0.012 | -0.004 | 0.027 | 1.000 |
No parece que haya mejorado mucho…Este estudio nos lleva a pensar que las variables continuas no van a aportar mucho a nuestro modelo de regresión, si bien puede suceder que en presencia de otras variables (por ejemplo las categóricas) la influencia de éstas pueda aumentar debido al distinto comportamiento (pendiente) que presenten con la respuesta en los distintos grupos que forma la variable categórica.
Por este motivo, es positivo complementar este análisis descriptivo con la información que dan los modelos sobre la importancia de variables, ya que dentro de modelo, la influencia de la variable se entiende a niveles constantes de todas las demás presentes en el mismo.
En cualquier caso, una vía que se puede explorar aquí para conseguir mayor influencia de las continuas es intentar tramificarlas (pasarlas a categóricas) de tal manera que se puedan descubrir patrones muy alejados d la linealidad que esta puedan contener. En el mundo de la tramificación existen dos estrategias fundamentales,
tramos ad-hoc para una mejor interpretabilidad (por ejemplo tramos de edad acordes a los que se utilizan en las estadísticas oficiales) en los que normalmente se realizan particiones más o menos equidistantes o basadas en un fundamento anterior.
tramos que maximicen cierta relación con la variable objetivo. Se pueden generar los puntos de corte de la variable a tramificar tales que los grupos creados sean lo más distintos entre sí en distribución de la variable respuesta, con lo que se “aseguran” tramos con capacidad de discriminar frente al objetivo. Suelen funcionar bien este tipo de tramificaciones pero tienen cierto peligro de sobreajustar a los datos de entrenamiento…Una de las formas de hacerlo es mediante árboles de regresión/clasificación.
Un árbol realiza particiones binarias de las variables con la premisa de encontrar las mayores diferencias entre los grupos que va formando, justo lo que buscamos!
# http://gnpalencia.org/optbinning/binning_continuous.html
from optbinning import ContinuousOptimalBinning
optb = ContinuousOptimalBinning(name='Azucar', dtype="numerical", max_n_bins=3)
optb.fit(imput_transf_cor['Azucar'].values, varObjCont)ContinuousOptimalBinning(max_n_bins=3, name='Azucar')In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
ContinuousOptimalBinning(max_n_bins=3, name='Azucar')
x_transform_bins = optb.transform(imput_transf_cor['Azucar'].values, metric="bins")
x_transform_bins## array(['(-inf, 44.35)', '(-inf, 44.35)', '(-inf, 44.35)', ...,
## '(-inf, 44.35)', '(-inf, 44.35)', '(-inf, 44.35)'], dtype=object)sns.violinplot(x=x_transform_bins,y=varObjCont,palette='viridis')
Podemos probar a realizar el binning de alguna (o todas) las variables continuas que no parezcan relevantes en el anterior análisis en pos de la mejor capacidad predictiva frente al modelo. Recomendable evaluar la salida sin restricciones del árbol y tal vez refinar con un máximo de categorías de salida en su caso. EL violin, revela las ligerísimas diferencias en distribución que se consiguen en este caso tramificando la variable Azucar.
En cualquier caso, es una opción interesante a considerar para la modelización predictiva.
Cuidado con pasarse de categorías!! Como ya sabemos las variables categóricas “consumen” grados de libertad y añaden k-1 parámetros al modelo, aumentando la complejidad del mismo y puediendo incurrir en problemas de robustez o generalización a nuevos datos. Por ello, hay que controlar el crecimiento de categorías de los factores.
Vamos a valorar la capacidad a priori de esta variable.
# ANOVA Azucar tramificada
vinosCompra['Azucar_rec'] = x_transform_bins
model = ols('Beneficio ~ C(Azucar_rec)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table## sum_sq df F PR(>F)
## C(Azucar_rec) 2.477983e+05 2.0 2.488539 0.083134
## Residual 2.486906e+08 4995.0 NaN NaNNo parece que haya una asociación relevante a nivel estadístico.
Son muchas las posibilidades sobre transformaciones de variables de cara a la modelización predictiva. Una estrategia interesante es aplicar múltiples vías de modelización con distintas posibilidades y comparar los resultados finales, decidiendo en base a la capacidad preditiva de la solución y su interpretabilidad.
Por el momento, nos centraremos en la creación de modelos de regresión para estimar el Beneficio de los vinos con Compra = 1 debido a esa mejor en la distribcuión de la respuesta.
Modelos de regresión lineal para la predicción del beneficio del vino
En esta sección se ajustan distintos modelos de regresión lineal para predecir el beneficio de los vinos con compra = 1. En primer lugar, tomamos la partición training (donde ajustamos el modelo) y test (donde probamos su capacidad).
Partición training-test
La filosofía de muchos de los modelos predictivos en python es la de tener por separado la matriz de predictores y el vector de variable objetivo ya que bajo este esquema están programados los procedimientos. Sin embargo, esto no es algo directo cuando existen variables categóricas en el dataset. En general, estas variables no son aceptadas tal cual en los modelos y hay que realizar un paso previo de creación de variables dummy generando lo que se conoce como la matriz de diseño del modelo.
Recordamos aquí lo que ya hemos comentado. Si tengo 4 categorías en un factor y pretendo generar dummies para todas, qué pasa?? Pues que tenemos una combinación lineal exacta de las variables en la matriz y por tanto no se puede invertir (algo necesario para la estimación por mínimos cuadrados de los parámetros del modelo) y todo se va al garete!! Por ello habría que tener en cuenta esto para generar solamente 3 dummies quedando “implícito” el efecto de la cuarta cuando todas ellas tomen el valor 0. Ok…y En este momento nos preguntamos, y donde va el efecto de esta variable? es decir, como se contabiliza en el modelo? Pues precisamente en nuestra querida constante \(\beta_0\) u ordenada en el origne de la ecuación de regresión. Vale, y donde está esta constante…pues depende..si generamos todo “a mano” es decir la extensión en k-1 dummies, tendríamos que añadir, así mismo esta constante (en realidad un vector de todo 1’s en la matriz de diseño).
# Función necesaria
from sklearn.model_selection import train_test_split
# Creamos 4 objetos: predictores para tr y tst y variable objetivo para tr y tst.
X_train, X_test, y_train, y_test = train_test_split(imputCompra, varObjCont, test_size=0.2, random_state=42)
# Comprobamos dimensiones
print('Training dataset shape:', X_train.shape, y_train.shape)## Training dataset shape: (3998, 16) (3998,)print('Testing dataset shape:', X_test.shape, y_test.shape)## Testing dataset shape: (1000, 16) (1000,)Vamos a explorar algunas posibilidades para automatizar estas cosas y evitar de momento la necesidad de generar explícitamente y mucho menos a mano, las matrices de diseño.
Algo muy guay que tiene R es la interfaz fórmula, mediante la cual podemos undicar el modelos que queramos ajustar de una forma muy visual. Por ejemplo, como hemos hecho en ANOVA, Beneficio en función de Región o Etiqueta o los efectos que queramos incluir…no tenemos más que poner ‘Beneficio ~ Efecto1 + Efecto2 +..+EfectoN’ y por debajo se estarán generando estas matrices de diseño con sus dummies en caso de necesidad y su constante y de todo.
Tenemos esta posibilidad en python para los modelos de regresión. Generemos primero el data_train con la variable objetivo dentro y posterioremente hacemos uso de la api para fórmulas de statmodels en su variante para regresión por mínimos cuadrados ordinarios.
# Genero el training con la objetivo dentro
data_train = X_train.join(y_train)# Importamos la api para fórmulas (en concreto ols para regresión)
from statsmodels.formula.api import ols
# Ajusto regresión de ejemplo
results = ols('Beneficio ~ Etiqueta + Acidez',data=data_train).fit()
results.summary()| Dep. Variable: | Beneficio | R-squared: | 0.370 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.370 |
| Method: | Least Squares | F-statistic: | 469.7 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:09:52 | Log-Likelihood: | -26401. |
| No. Observations: | 3998 | AIC: | 5.281e+04 |
| Df Residuals: | 3992 | BIC: | 5.285e+04 |
| Df Model: | 5 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 288.7290 | 13.962 | 20.680 | 0.000 | 261.356 | 316.102 |
| Etiqueta[T.M] | 138.7840 | 14.979 | 9.265 | 0.000 | 109.416 | 168.152 |
| Etiqueta[T.R] | 291.5636 | 14.522 | 20.077 | 0.000 | 263.092 | 320.036 |
| Etiqueta[T.B] | 452.1047 | 15.051 | 30.037 | 0.000 | 422.596 | 481.614 |
| Etiqueta[T.MB] | 585.2175 | 20.101 | 29.115 | 0.000 | 545.809 | 624.626 |
| Acidez | -7.0863 | 3.685 | -1.923 | 0.055 | -14.311 | 0.139 |
| Omnibus: | 131.498 | Durbin-Watson: | 1.991 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 143.931 |
| Skew: | 0.461 | Prob(JB): | 5.57e-32 |
| Kurtosis: | 3.122 | Cond. No. | 13.8 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Se peude comprobar que el funcionamiento es adecuado en tanto en cuanto la variable categórica Etiqueta aparece con 5-1 = 4 efectos en el modelo (siendo el restante ‘MM’ recogido por la constante \(\beta_0\)). Además, vemos que se conserva el órden de las categorías que hemos especificado y también que el efetco tiene pinta de monótono creciente con etiqueta puesto que los betas van aumentando con el aumento de etiqueta. Si nos ponemos a valorar, pues es un modelo bastante malo con un R2 de 0,37 pero todos los parámetros son significativos a nivel estadístico (es decir, su parámetro es distinto de 0 o el IC no contiene al valor 0). Ajuste muy pobre con significación paramétrica.
Lo que siempre haremos es un modelo completo de referencia para valorar la capacidad (sobre training) de este y tomarla como base para la generación de modelos más sencillos que mantengan valor predictivo en training y lo mejoren en test.
Modelo completo de referencia
Para utilizar la interfaz fórmula proporcionada por la api, tenemos que especificar todos los efectos…cosa que no me gusta nada (R tiene la posibilidad de poner ‘Beneficio ~ .’ y se sobreentiende que quieres todos los efectos) pero la vida es durilla a veces jeje.
Venga pues vamos a facilitarnos la vida un poco y generamos una función que concatene todos los efectos como string e incluso podamos eliminar algunos a placer, obteniendo la fórmula de manera automática.
# Función para generar la fórmula por larga que sea
def ols_formula(df, dependent_var, *excluded_cols):
df_columns = list(df.columns.values)
df_columns.remove(dependent_var)
for col in excluded_cols:
df_columns.remove(col)
return dependent_var + ' ~ ' + ' + '.join(df_columns)
# Aplicamos a fórmula de modelo completo
form=ols_formula(data_train,'Beneficio')
form## 'Beneficio ~ Acidez + AcidoCitrico + Azucar + CloruroSodico + Densidad + pH + Sulfatos + Alcohol + CalifProductor + PrecioBotella + Etiqueta + Clasificacion + Region + prop_missings + aleatorio + aleatorio2'# Ajusto regresión según fórmula completa
modeloC = ols(form,data=data_train).fit()
modeloC.summary()| Dep. Variable: | Beneficio | R-squared: | 0.435 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.431 |
| Method: | Least Squares | F-statistic: | 122.2 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:10:01 | Log-Likelihood: | -26186. |
| No. Observations: | 3998 | AIC: | 5.242e+04 |
| Df Residuals: | 3972 | BIC: | 5.259e+04 |
| Df Model: | 25 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 273.6801 | 107.982 | 2.534 | 0.011 | 61.974 | 485.386 |
| Etiqueta[T.M] | 123.1336 | 14.323 | 8.597 | 0.000 | 95.052 | 151.215 |
| Etiqueta[T.R] | 252.3407 | 14.069 | 17.937 | 0.000 | 224.758 | 279.923 |
| Etiqueta[T.B] | 383.3082 | 14.824 | 25.857 | 0.000 | 354.244 | 412.372 |
| Etiqueta[T.MB] | 495.2993 | 19.721 | 25.115 | 0.000 | 456.635 | 533.963 |
| Clasificacion[T.*] | 17.9532 | 9.192 | 1.953 | 0.051 | -0.068 | 35.974 |
| Clasificacion[T.**] | 63.0280 | 8.857 | 7.117 | 0.000 | 45.664 | 80.392 |
| Clasificacion[T.***] | 127.4193 | 9.827 | 12.966 | 0.000 | 108.152 | 146.686 |
| Clasificacion[T.****] | 205.0624 | 13.503 | 15.187 | 0.000 | 178.590 | 231.535 |
| Region[T.2.0] | 12.1224 | 6.623 | 1.830 | 0.067 | -0.863 | 25.107 |
| Region[T.3.0] | 0.9415 | 6.565 | 0.143 | 0.886 | -11.929 | 13.812 |
| prop_missings[T.7.6923076923076925] | 3.0431 | 6.504 | 0.468 | 0.640 | -9.708 | 15.794 |
| prop_missings[T.16.666666666666664] | 21.9144 | 17.222 | 1.272 | 0.203 | -11.850 | 55.679 |
| prop_missings[T.27.27272727272727] | -49.3978 | 64.388 | -0.767 | 0.443 | -175.634 | 76.838 |
| Acidez | -7.2167 | 3.510 | -2.056 | 0.040 | -14.099 | -0.334 |
| AcidoCitrico | -1.6689 | 3.222 | -0.518 | 0.605 | -7.986 | 4.648 |
| Azucar | -0.1400 | 0.079 | -1.778 | 0.076 | -0.294 | 0.014 |
| CloruroSodico | -5.9563 | 8.484 | -0.702 | 0.483 | -22.590 | 10.678 |
| Densidad | -61.1204 | 106.290 | -0.575 | 0.565 | -269.509 | 147.269 |
| pH | 3.8509 | 4.030 | 0.956 | 0.339 | -4.050 | 11.752 |
| Sulfatos | 0.7644 | 2.885 | 0.265 | 0.791 | -4.892 | 6.421 |
| Alcohol | 4.3407 | 0.751 | 5.779 | 0.000 | 2.868 | 5.813 |
| CalifProductor | -6.3493 | 2.502 | -2.538 | 0.011 | -11.254 | -1.445 |
| PrecioBotella | -0.3459 | 1.844 | -0.188 | 0.851 | -3.961 | 3.270 |
| aleatorio | -3.4986 | 9.380 | -0.373 | 0.709 | -21.888 | 14.891 |
| aleatorio2 | 12.1311 | 9.311 | 1.303 | 0.193 | -6.123 | 30.386 |
| Omnibus: | 113.210 | Durbin-Watson: | 1.965 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 122.531 |
| Skew: | 0.429 | Prob(JB): | 2.47e-27 |
| Kurtosis: | 3.029 | Cond. No. | 1.96e+03 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.96e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Modelo completo bastante sobreparametrizado por lo que se ve, con muchos efectos no significativos cuya aportación al R2 será seguramente despreciable y puedan ser eliminados y con una capacidad bastante baja con R2 de 0.43…mal augurio. Conjunto de datos dificil de modelizar.. Pero hay que tirar palante con lo que tenemos y luchar!
Importacia de las variables
Vamos a tratar de extraer la importancia de las variables según su aportación al R2 total del modelo. Con la función relativeImp podemos realizar esta labor. Sin embargo, de nuevo u pequeño obstáculo…solamente funciona con variables continuas..vaya!
Probamos.
from relativeImp import relativeImp
# Nombres de las variables continuas
names=X_train.select_dtypes(include=np.number).columns.tolist()
# Calcular importancia relativa de las variable continuas en modelo (solo con variables continuas!)
df_results = relativeImp(data_train, outcomeName = 'Beneficio', driverNames = names)
df_results.sort_values(by='normRelaImpt', ascending=False)## driver rawRelaImpt normRelaImpt
## 7 Alcohol 0.006235 52.287319
## 8 CalifProductor 0.002157 18.086350
## 0 Acidez 0.001874 15.715050
## 11 aleatorio2 0.000492 4.129008
## 2 Azucar 0.000465 3.898723
## 5 pH 0.000180 1.507559
## 9 PrecioBotella 0.000158 1.328349
## 4 Densidad 0.000151 1.266428
## 6 Sulfatos 0.000074 0.623407
## 3 CloruroSodico 0.000065 0.543347
## 10 aleatorio 0.000037 0.313968
## 1 AcidoCitrico 0.000036 0.300494Con esta funcionalidad siempre podemos calcular la importacia relativa según la aportación al R2 global de cada una de las variables de un modelo especifico. Sin embargo, la función no acepta variables categóricas por lo que, para poder tener información sobre la importacia del modelo completo, es necesario realizar la extensión en dummies (algo que sucede en el interior del modelo de regresión con la api de fórmula).
En el modelo solamente con variables continuas (con un R2 lamentable como de 0.01) las más importantes son Alcohol, CalifProductor y Acidez.
Lo que nos gustaría es poder ver esto para el modelo completo incluyendo las nominales. Para ello, no nos queda otra alternativa que generar de manera explícita las matrices de diseño del modelo.
El paquete pasty resulta muy útil para generar las matrices de diseño de un modelo basadas en una fórmula concreta. Esto nos permitirá utilizar cualquier modelo de predicción de python que no acepte directamente la interfaz fórmula. Es evidente que podríamos hacer esta operación “a mano”, generando la extensión en dummies de los factores y la posterior adición de una constante \(\beta_0\) para que se refleje en el modelo, pero pienso que es más intuitivo utilizar la fórmula y así evitar errores.
import patsy
# Generamos las matrices de diseño según la fórmula de modelo completo
y, X = patsy.dmatrices(form, data_train, return_type='dataframe')
# Ahora podemos aplicar la función "oficial" de statmodels OLS (con formato y,X)
model=sm.OLS(y,X).fit()
model.summary()| Dep. Variable: | Beneficio | R-squared: | 0.435 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.431 |
| Method: | Least Squares | F-statistic: | 122.2 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:10:09 | Log-Likelihood: | -26186. |
| No. Observations: | 3998 | AIC: | 5.242e+04 |
| Df Residuals: | 3972 | BIC: | 5.259e+04 |
| Df Model: | 25 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 273.6801 | 107.982 | 2.534 | 0.011 | 61.974 | 485.386 |
| Etiqueta[T.M] | 123.1336 | 14.323 | 8.597 | 0.000 | 95.052 | 151.215 |
| Etiqueta[T.R] | 252.3407 | 14.069 | 17.937 | 0.000 | 224.758 | 279.923 |
| Etiqueta[T.B] | 383.3082 | 14.824 | 25.857 | 0.000 | 354.244 | 412.372 |
| Etiqueta[T.MB] | 495.2993 | 19.721 | 25.115 | 0.000 | 456.635 | 533.963 |
| Clasificacion[T.*] | 17.9532 | 9.192 | 1.953 | 0.051 | -0.068 | 35.974 |
| Clasificacion[T.**] | 63.0280 | 8.857 | 7.117 | 0.000 | 45.664 | 80.392 |
| Clasificacion[T.***] | 127.4193 | 9.827 | 12.966 | 0.000 | 108.152 | 146.686 |
| Clasificacion[T.****] | 205.0624 | 13.503 | 15.187 | 0.000 | 178.590 | 231.535 |
| Region[T.2.0] | 12.1224 | 6.623 | 1.830 | 0.067 | -0.863 | 25.107 |
| Region[T.3.0] | 0.9415 | 6.565 | 0.143 | 0.886 | -11.929 | 13.812 |
| prop_missings[T.7.6923076923076925] | 3.0431 | 6.504 | 0.468 | 0.640 | -9.708 | 15.794 |
| prop_missings[T.16.666666666666664] | 21.9144 | 17.222 | 1.272 | 0.203 | -11.850 | 55.679 |
| prop_missings[T.27.27272727272727] | -49.3978 | 64.388 | -0.767 | 0.443 | -175.634 | 76.838 |
| Acidez | -7.2167 | 3.510 | -2.056 | 0.040 | -14.099 | -0.334 |
| AcidoCitrico | -1.6689 | 3.222 | -0.518 | 0.605 | -7.986 | 4.648 |
| Azucar | -0.1400 | 0.079 | -1.778 | 0.076 | -0.294 | 0.014 |
| CloruroSodico | -5.9563 | 8.484 | -0.702 | 0.483 | -22.590 | 10.678 |
| Densidad | -61.1204 | 106.290 | -0.575 | 0.565 | -269.509 | 147.269 |
| pH | 3.8509 | 4.030 | 0.956 | 0.339 | -4.050 | 11.752 |
| Sulfatos | 0.7644 | 2.885 | 0.265 | 0.791 | -4.892 | 6.421 |
| Alcohol | 4.3407 | 0.751 | 5.779 | 0.000 | 2.868 | 5.813 |
| CalifProductor | -6.3493 | 2.502 | -2.538 | 0.011 | -11.254 | -1.445 |
| PrecioBotella | -0.3459 | 1.844 | -0.188 | 0.851 | -3.961 | 3.270 |
| aleatorio | -3.4986 | 9.380 | -0.373 | 0.709 | -21.888 | 14.891 |
| aleatorio2 | 12.1311 | 9.311 | 1.303 | 0.193 | -6.123 | 30.386 |
| Omnibus: | 113.210 | Durbin-Watson: | 1.965 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 122.531 |
| Skew: | 0.429 | Prob(JB): | 2.47e-27 |
| Kurtosis: | 3.029 | Cond. No. | 1.96e+03 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.96e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
# Ajusto regresión de ejemplo
results = ols('Beneficio ~ aleatorio2 + Etiqueta',data=data_train).fit()
results.summary()| Dep. Variable: | Beneficio | R-squared: | 0.370 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.369 |
| Method: | Least Squares | F-statistic: | 469.3 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:10:14 | Log-Likelihood: | -26402. |
| No. Observations: | 3998 | AIC: | 5.282e+04 |
| Df Residuals: | 3992 | BIC: | 5.285e+04 |
| Df Model: | 5 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 278.1061 | 14.725 | 18.887 | 0.000 | 249.237 | 306.976 |
| Etiqueta[T.M] | 139.7310 | 14.976 | 9.330 | 0.000 | 110.369 | 169.093 |
| Etiqueta[T.R] | 292.7666 | 14.514 | 20.172 | 0.000 | 264.312 | 321.222 |
| Etiqueta[T.B] | 453.0996 | 15.044 | 30.118 | 0.000 | 423.605 | 482.595 |
| Etiqueta[T.MB] | 586.8203 | 20.095 | 29.202 | 0.000 | 547.422 | 626.218 |
| aleatorio2 | 14.7823 | 9.784 | 1.511 | 0.131 | -4.401 | 33.965 |
| Omnibus: | 131.244 | Durbin-Watson: | 1.993 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 143.626 |
| Skew: | 0.460 | Prob(JB): | 6.49e-32 |
| Kurtosis: | 3.122 | Cond. No. | 14.3 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Podemos comprobar que el modelo es exactamente el mismo que el ajustado mediante la fórmula y la función ols, pero ahora resulta que las matrices de diseño se especifican explicitamente con lo que podemos utilizar la importancia relativa de variables ya que no existen factores. Veamos..
# Nombres de predictores (en modo dummy) donde quitamos la constante
names=X.columns.tolist()[1:]
# Calculamos importancia relativa
df_results = relativeImp(X.join(y), outcomeName = 'Beneficio', driverNames = names)
# Ordenamos valores
df_results.sort_values(by='normRelaImpt', ascending=False)## driver rawRelaImpt normRelaImpt
## 2 Etiqueta[T.B] 0.141921 32.645204
## 3 Etiqueta[T.MB] 0.072831 16.752895
## 7 Clasificacion[T.****] 0.057918 13.322638
## 0 Etiqueta[T.M] 0.048181 11.082912
## 6 Clasificacion[T.***] 0.046103 10.604822
## 1 Etiqueta[T.R] 0.034920 8.032400
## 4 Clasificacion[T.*] 0.014820 3.408889
## 5 Clasificacion[T.**] 0.006730 1.548010
## 20 Alcohol 0.005406 1.243454
## 21 CalifProductor 0.001583 0.364180
## 13 Acidez 0.001324 0.304557
## 8 Region[T.2.0] 0.000728 0.167488
## 15 Azucar 0.000522 0.119964
## 11 prop_missings[T.16.666666666666664] 0.000413 0.095110
## 24 aleatorio2 0.000349 0.080353
## 12 prop_missings[T.27.27272727272727] 0.000261 0.059923
## 18 pH 0.000213 0.048928
## 9 Region[T.3.0] 0.000111 0.025580
## 16 CloruroSodico 0.000104 0.023823
## 17 Densidad 0.000099 0.022765
## 14 AcidoCitrico 0.000056 0.012985
## 22 PrecioBotella 0.000046 0.010545
## 19 Sulfatos 0.000045 0.010264
## 23 aleatorio 0.000034 0.007854
## 10 prop_missings[T.7.6923076923076925] 0.000019 0.004455# Gráfico de importancia relativa en base al R2
px.bar(df_results,x='normRelaImpt',y='driver',title='Importancia relativa por aportación al R2').update_yaxes(categoryorder="total ascending").show()Ahora si que si! Obtenemos un gráfico de importancia de las variables del modelo completo con el efecto de las categorías de las nominales por separado. Claramente lo que se podía esperar, las categorías de Clasificación y Etiqueta son las de mayor importancia para el modelo.
Vamos a generar un proceso backward eliminando variables secuencialmente según el p-valor.
# Proceso backward de eliminación de efectos según p-valor
form2=ols_formula(data_train,'Beneficio','prop_missings','aleatorio','aleatorio2',
'PrecioBotella','Sulfatos','AcidoCitrico','Densidad',
'CloruroSodico','pH','Region','Azucar')
# Ajusto regresión sin prop_missings
modeloC2 = ols(form2,data=data_train).fit()
modeloC2.summary()| Dep. Variable: | Beneficio | R-squared: | 0.433 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.431 |
| Method: | Least Squares | F-statistic: | 276.5 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:10:25 | Log-Likelihood: | -26193. |
| No. Observations: | 3998 | AIC: | 5.241e+04 |
| Df Residuals: | 3986 | BIC: | 5.248e+04 |
| Df Model: | 11 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 232.5770 | 17.809 | 13.059 | 0.000 | 197.661 | 267.493 |
| Etiqueta[T.M] | 123.9191 | 14.299 | 8.666 | 0.000 | 95.885 | 151.954 |
| Etiqueta[T.R] | 253.0856 | 14.042 | 18.023 | 0.000 | 225.555 | 280.616 |
| Etiqueta[T.B] | 384.7283 | 14.798 | 26.000 | 0.000 | 355.717 | 413.740 |
| Etiqueta[T.MB] | 496.4143 | 19.690 | 25.211 | 0.000 | 457.810 | 535.019 |
| Clasificacion[T.*] | 17.7085 | 9.173 | 1.930 | 0.054 | -0.276 | 35.693 |
| Clasificacion[T.**] | 62.7754 | 8.842 | 7.100 | 0.000 | 45.440 | 80.111 |
| Clasificacion[T.***] | 126.8453 | 9.814 | 12.925 | 0.000 | 107.604 | 146.087 |
| Clasificacion[T.****] | 204.3448 | 13.491 | 15.147 | 0.000 | 177.895 | 230.794 |
| Acidez | -6.9836 | 3.501 | -1.995 | 0.046 | -13.848 | -0.119 |
| Alcohol | 4.4198 | 0.750 | 5.895 | 0.000 | 2.950 | 5.890 |
| CalifProductor | -6.8412 | 2.481 | -2.757 | 0.006 | -11.705 | -1.977 |
| Omnibus: | 114.292 | Durbin-Watson: | 1.967 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 123.837 |
| Skew: | 0.431 | Prob(JB): | 1.29e-27 |
| Kurtosis: | 3.016 | Cond. No. | 135. |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
from sklearn.metrics import mean_squared_error, r2_score
# Predicciones para test modelo completo
vinos_y_pred = modeloC.predict(X_test)
# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, vinos_y_pred)))
# The coefficient of determination: 1 is perfect prediction## Mean squared error: 161.55print("Coefficient of determination: %.2f" % r2_score(y_test, vinos_y_pred))
# Gráfico de real frente a prediccion para alguna variable
#plt.scatter(X_test.Acidez, y_test, color="black")
#plt.scatter(X_test.Acidez, vinos_y_pred, color="blue")## Coefficient of determination: 0.44from sklearn.metrics import mean_squared_error, r2_score
# Predicciones para test modelo C2
vinos_y_pred = modeloC2.predict(X_test)
# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, vinos_y_pred)))
# The coefficient of determination: 1 is perfect prediction## Mean squared error: 161.31print("Coefficient of determination: %.2f" % r2_score(y_test, vinos_y_pred))
# Gráfico de real frente a prediccion para alguna variable
#plt.scatter(X_test.Acidez, y_test, color="black")
#plt.scatter(X_test.Acidez, vinos_y_pred, color="blue")## Coefficient of determination: 0.44Ajuste por validación cruzada repetida
En esta sección vamos a valorar el ajuste de modelos por el proceso de validación cruzada repetida, en el que se generarán n_splits particiones del archivo, repitiendo el proceso con distintas semillas n_repeats veces. De esta forma, obtenemos n_splits x n_repeats modelos por fórmula especificada, pudiendo así promediar los resultados obtenidos bajo muchas muestras de training-test.
En este esquema, todas las observaciones del archivo son utilizadas en algún training para el ajuste del modelo yen algún test para su evaluación pero nunca simultáneamente. Así el modelo tiene instancias cambiantes para el ajuste de parámetros y para la predicción.
# from sklearn.model_selection import cross_val_score
# from sklearn.model_selection import RepeatedKFold
# from sklearn.linear_model import LinearRegression
model = LinearRegression()
#model = sm.OLS(y,X)
# Establecemos esquema de validación fijando random_state (reproducibilidad)
cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=12345)
# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X, y, cv=cv)
# Sesgo y varianza
print('Coeficioente de detrminación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))## Coeficioente de detrminación R2: 0.425 (0.023)sns.violinplot(y=scores,palette='viridis')
Modelo final backward
# Generamos las matrices de diseño según la fórmula de modelo completo
y2, X2 = patsy.dmatrices(form2, vinosCompra, return_type='dataframe')
# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X2, y2, cv=cv)
# Sesgo y varianza
print('Coeficioente de detrminación R2 Modelo Final Backward: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))## Coeficioente de detrminación R2 Modelo Final Backward: 0.430 (0.020)sns.violinplot(y=scores,palette='viridis')
Modelo con interacción Etiqueta-Clasificación
Siendo variables que miden de alguna forma la calidad del vino, es posible que exista cierta interacción entre ellas, es decir que la pertenencia conjunta a determinada etiqueta y clasificación tenga efectos diferenciadores en el beneficio que produce el vino.
form3 = 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta*Clasificacion'
# Ajusto regresión sin prop_missings
modeloC3 = ols(form3,data=data_train).fit()
modeloC3.summary()| Dep. Variable: | Beneficio | R-squared: | 0.439 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.435 |
| Method: | Least Squares | F-statistic: | 119.6 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:10:51 | Log-Likelihood: | -26170. |
| No. Observations: | 3998 | AIC: | 5.239e+04 |
| Df Residuals: | 3971 | BIC: | 5.256e+04 |
| Df Model: | 26 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 214.0458 | 25.322 | 8.453 | 0.000 | 164.401 | 263.690 |
| Etiqueta[T.M] | 123.3191 | 25.968 | 4.749 | 0.000 | 72.408 | 174.231 |
| Etiqueta[T.R] | 265.0898 | 26.056 | 10.174 | 0.000 | 214.006 | 316.174 |
| Etiqueta[T.B] | 481.7667 | 33.008 | 14.596 | 0.000 | 417.053 | 546.481 |
| Etiqueta[T.MB] | 694.7219 | 60.959 | 11.396 | 0.000 | 575.207 | 814.237 |
| Clasificacion[T.*] | 48.5567 | 30.214 | 1.607 | 0.108 | -10.680 | 107.793 |
| Clasificacion[T.**] | 91.3712 | 40.012 | 2.284 | 0.022 | 12.926 | 169.816 |
| Clasificacion[T.***] | 95.6293 | 60.987 | 1.568 | 0.117 | -23.940 | 215.198 |
| Clasificacion[T.****] | 124.0220 | 17.909 | 6.925 | 0.000 | 88.911 | 159.133 |
| Etiqueta[T.M]:Clasificacion[T.*] | -8.2539 | 33.620 | -0.246 | 0.806 | -74.168 | 57.660 |
| Etiqueta[T.R]:Clasificacion[T.*] | -22.9217 | 33.481 | -0.685 | 0.494 | -88.563 | 42.720 |
| Etiqueta[T.B]:Clasificacion[T.*] | -144.7659 | 40.925 | -3.537 | 0.000 | -225.003 | -64.529 |
| Etiqueta[T.MB]:Clasificacion[T.*] | -112.5326 | 87.554 | -1.285 | 0.199 | -284.189 | 59.123 |
| Etiqueta[T.M]:Clasificacion[T.**] | -13.2018 | 42.632 | -0.310 | 0.757 | -96.785 | 70.381 |
| Etiqueta[T.R]:Clasificacion[T.**] | -19.0431 | 42.242 | -0.451 | 0.652 | -101.861 | 63.774 |
| Etiqueta[T.B]:Clasificacion[T.**] | -113.6291 | 47.330 | -2.401 | 0.016 | -206.423 | -20.835 |
| Etiqueta[T.MB]:Clasificacion[T.**] | -224.5732 | 74.642 | -3.009 | 0.003 | -370.914 | -78.233 |
| Etiqueta[T.M]:Clasificacion[T.***] | 74.5139 | 64.540 | 1.155 | 0.248 | -52.021 | 201.049 |
| Etiqueta[T.R]:Clasificacion[T.***] | 26.5701 | 62.673 | 0.424 | 0.672 | -96.304 | 149.444 |
| Etiqueta[T.B]:Clasificacion[T.***] | -38.5356 | 66.105 | -0.583 | 0.560 | -168.139 | 91.068 |
| Etiqueta[T.MB]:Clasificacion[T.***] | -161.0492 | 85.715 | -1.879 | 0.060 | -329.099 | 7.000 |
| Etiqueta[T.M]:Clasificacion[T.****] | 160.8211 | 44.287 | 3.631 | 0.000 | 73.994 | 247.648 |
| Etiqueta[T.R]:Clasificacion[T.****] | 101.9520 | 25.367 | 4.019 | 0.000 | 52.218 | 151.686 |
| Etiqueta[T.B]:Clasificacion[T.****] | -4.1577 | 27.740 | -0.150 | 0.881 | -58.544 | 50.229 |
| Etiqueta[T.MB]:Clasificacion[T.****] | -134.5934 | 52.123 | -2.582 | 0.010 | -236.784 | -32.403 |
| Acidez | -7.2240 | 3.500 | -2.064 | 0.039 | -14.086 | -0.362 |
| Alcohol | 4.4645 | 0.748 | 5.967 | 0.000 | 2.998 | 5.931 |
| CalifProductor | -6.5027 | 2.477 | -2.626 | 0.009 | -11.358 | -1.647 |
| Omnibus: | 108.268 | Durbin-Watson: | 1.966 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 116.838 |
| Skew: | 0.419 | Prob(JB): | 4.26e-26 |
| Kurtosis: | 2.992 | Cond. No. | 1.05e+16 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 4.82e-27. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
# Generamos las matrices de diseño según la fórmula de modelo completo
y3, X3 = patsy.dmatrices(form3, vinosCompra, return_type='dataframe')
# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X3, y3, cv=cv)
# Sesgo y varianza
print('Coeficioente de detrminación R2 Modelo Final Backward: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))## Coeficioente de detrminación R2 Modelo Final Backward: 0.422 (0.027)sns.violinplot(y=scores,palette='viridis')
Llegados a este punto se observa que, si bien el modelo con la interacción es bastante bueno, tiene un cruce de categorías de la interacción que no contiene observaciones en data_train para poder ser estimado su parámetro…esto es un fallo en la especificación del modelo y será difícil mantenerlo tal cual puesto que a nivel interpretativo carece de la robustez necesaria.
Una estrategia que se puede seguir es unir alguna de las categorías implicadas en esta interacción peligrosa. En este caso podemos unir las categorías de clasificación más elevadas (se puede probar trabajando sobre las etiquetas!!).
Cuando hacemos lago así, hay que tener la precaución de hacer el cambio en el dataset completo y volver a generar la partición para que los datos se actualicen en training y también en test.
Actualizaremos el último modelo con esta unión de categorías.
imputCompra['Etiqueta_r'] = vinosCompra.Etiqueta_r
imputCompra['Azucar_rec'] = vinosCompra.Azucar_rec
# Creamos 4 objetos: predictores para tr y tst y variable objetivo para tr y tst.
X_train, X_test, y_train, y_test = train_test_split(imputCompra, varObjCont, test_size=0.2, random_state=42)
# Genero el training con la objetivo dentro
data_train_r = X_train.join(y_train)
form4 = 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta_r*Clasificacion'
# Generamos las matrices de diseño según la fórmula de modelo completo
y4, X4 = patsy.dmatrices(form4, vinosCompra, return_type='dataframe')
# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X4, y4, cv=cv)
# Sesgo y varianza
print('Coeficioente de detrminación R2 Modelo Final Backward interacción Rec: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))## Coeficioente de detrminación R2 Modelo Final Backward interacción Rec: 0.421 (0.019)sns.violinplot(y=scores,palette='viridis')
Comparación por validación cruzada
Una vez explorados varios modelos manuales, es bueno comprobar sus capacidades en un esquema de validación cruzada repetida con la intención de comprobar su estabilidad ante el remuestreo.
El objetivo fundamental en esta parte es seleccionar el mejor modelo en relación sesgo-varianza y complejidad.
Vamos a generar una función que automatice un poco el proceso de comparación por validación cruzada, de tal forma que pueda aplicar sobre cualquier fórmula epecificada. Luego será cuestión de enlistar todas las fómulas que queramos y pasar la función con un map().
# Función para comparación por validación cruzada
def cross_val(formula, data, seed=12345):
# Generamos las matrices de diseño según la fórmula de modelo completo
y, X = patsy.dmatrices(formula, data, return_type='dataframe')
# Establecemos esquema de validación fijando random_state (reproducibilidad)
cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=seed)
# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X, y, cv=cv)
# Sesgo y varianza
print('Modelo: ' + formula)
print('Coeficiente de determinación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
#sns.violinplot(y=scores,palette='viridis')
return(scores)form5 = 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion + Azucar_rec'
form6 = 'Beneficio ~ Alcohol + CalifProductor + Etiqueta + Clasificacion'
form7 = 'Beneficio ~ CalifProductor + Etiqueta + Clasificacion + Azucar_rec'
form8 = 'Beneficio ~ Etiqueta + Clasificacion + Azucar_rec'
# Creamos lista de fórmulas
list_form = [form,form2,form3,form4,form5,form6,form7,form8]
list_form
# Aplicamos a toda la lista la función creada (devuelve un dataframe pero está transpuesto)## ['Beneficio ~ Acidez + AcidoCitrico + Azucar + CloruroSodico + Densidad + pH + Sulfatos + Alcohol + CalifProductor + PrecioBotella + Etiqueta + Clasificacion + Region + prop_missings + aleatorio + aleatorio2', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta*Clasificacion', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta_r*Clasificacion', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion + Azucar_rec', 'Beneficio ~ Alcohol + CalifProductor + Etiqueta + Clasificacion', 'Beneficio ~ CalifProductor + Etiqueta + Clasificacion + Azucar_rec', 'Beneficio ~ Etiqueta + Clasificacion + Azucar_rec']list_res = pd.DataFrame(map(lambda x: cross_val(x,vinosCompra, seed=2022),list_form))
# Trasnponer dataframe y pasar de wide a long (creando un factor variable con el nombre de cada fórmula de la lista[0,1,2,3])## Modelo: Beneficio ~ Acidez + AcidoCitrico + Azucar + CloruroSodico + Densidad + pH + Sulfatos + Alcohol + CalifProductor + PrecioBotella + Etiqueta + Clasificacion + Region + prop_missings + aleatorio + aleatorio2
## Coeficiente de determinación R2: 0.429 (0.019)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion
## Coeficiente de determinación R2: 0.431 (0.018)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta*Clasificacion
## Coeficiente de determinación R2: 0.426 (0.018)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta_r*Clasificacion
## Coeficiente de determinación R2: 0.421 (0.019)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion + Azucar_rec
## Coeficiente de determinación R2: 0.431 (0.018)
## Modelo: Beneficio ~ Alcohol + CalifProductor + Etiqueta + Clasificacion
## Coeficiente de determinación R2: 0.430 (0.018)
## Modelo: Beneficio ~ CalifProductor + Etiqueta + Clasificacion + Azucar_rec
## Coeficiente de determinación R2: 0.426 (0.018)
## Modelo: Beneficio ~ Etiqueta + Clasificacion + Azucar_rec
## Coeficiente de determinación R2: 0.425 (0.019)results = list_res.T.melt()
results.columns = ['Modelo','R2']
results.head()## Modelo R2
## 0 0 0.394846
## 1 0 0.453454
## 2 0 0.388224
## 3 0 0.427713
## 4 0 0.457318Una vez generados los 100 modelos por fórmula con distintas particiones de trainig test bajo un esquema de validación cruzada repetida 5 Folds 20 Repeats, juntamos los resultados de los “resamples” en un dataset para generar el gráfico de boxplots de las distribuciones de R2 a través de los modelos y valorar la relación sesgo-varianza.
# Boxplot paralelo para comparar
sns.boxplot(x='Modelo',y='R2',data=results,palette='viridis')
En este punto hay que decidir un modelo final.
Todo indica que el mejor modelo en todos los sentidos es el modelo final backwar sin interacción. Modelo sencillo y con capacidad predictiva mayor en términos generales que sus competidores que además son más complejos. Claramente nos quedamos con el indicado como 1 que es el modelo de fórmula form2.
Ajuste del modelo final. Interpretación de parámetros
Una vez escogido el modelo final. Hay que interpretar sus parámetros teniendo en cuenta que la estimación más robusta de los mismos será la que se obtiene utilizando el datset completo y no solamente las instancias del conjunto de training.
De esta forma, se ajusta el modelo final a los datos completos.
# Ajusto regresión sin prop_missings
modelo_final = ols(form2,data=vinosCompra).fit()
modelo_final.summary()| Dep. Variable: | Beneficio | R-squared: | 0.435 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.433 |
| Method: | Least Squares | F-statistic: | 348.5 |
| Date: | do., 30 oct. 2022 | Prob (F-statistic): | 0.00 |
| Time: | 21:11:23 | Log-Likelihood: | -32696. |
| No. Observations: | 4998 | AIC: | 6.542e+04 |
| Df Residuals: | 4986 | BIC: | 6.549e+04 |
| Df Model: | 11 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 237.4493 | 15.767 | 15.060 | 0.000 | 206.539 | 268.359 |
| Etiqueta[T.M] | 123.4425 | 12.656 | 9.753 | 0.000 | 98.631 | 148.254 |
| Etiqueta[T.R] | 255.8916 | 12.431 | 20.585 | 0.000 | 231.522 | 280.261 |
| Etiqueta[T.B] | 385.2418 | 13.117 | 29.370 | 0.000 | 359.527 | 410.957 |
| Etiqueta[T.MB] | 494.7721 | 17.662 | 28.014 | 0.000 | 460.147 | 529.397 |
| Clasificacion[T.*] | 20.7499 | 8.105 | 2.560 | 0.010 | 4.860 | 36.640 |
| Clasificacion[T.**] | 66.6822 | 7.829 | 8.518 | 0.000 | 51.334 | 82.030 |
| Clasificacion[T.***] | 124.9397 | 8.658 | 14.431 | 0.000 | 107.967 | 141.913 |
| Clasificacion[T.****] | 208.6767 | 11.971 | 17.433 | 0.000 | 185.209 | 232.144 |
| Acidez | -7.3930 | 3.107 | -2.379 | 0.017 | -13.485 | -1.301 |
| Alcohol | 4.1374 | 0.665 | 6.218 | 0.000 | 2.833 | 5.442 |
| CalifProductor | -8.4882 | 2.216 | -3.831 | 0.000 | -12.832 | -4.145 |
| Omnibus: | 133.624 | Durbin-Watson: | 2.006 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 144.045 |
| Skew: | 0.415 | Prob(JB): | 5.26e-32 |
| Kurtosis: | 3.045 | Cond. No. | 136. |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Conclusiones que podemos sacar para aquellos vinos con Compra = 1, es decir beneficio distinto de 0:
Los vinos con Clasificación ** ** tienen un beneficio estimado medio 208.67 unidades superior a los vinos de clasificación Desconocida. El rango de variación esperado es [185.209, 232.144].
El beneficio estimado medio de los vinos de Clasificación * es 20.74 unidades superior al de los vinos con Clasificación Deconocida. El rango de variación esperado es [4.860, 36.640].
El beneficio estimado medio de los vinos de Etiqueta MB, es 494.77 unidades mayor que el de los vinos con etiqueta MM. El rango de variación esperado es [460.147, 529.397].
El aumento unitario de la Calificación del Productor, produce una disminución del beneficio estimado del vino de 8.488 unidades. Puede resutar sorprendente pero hay que tener en cuenta que la variación es a un nivel marginal, con lo cual se valora qué pasa con el aumento unitario de la calificación del propio productor para vinos con la misma clasificación, etiqueta y niveles de alcohol y acidez.
El aumento unitario en el nivel de Acidez produce una *disminución del Beneficio esperado de 7.39 unidades**
El aumento unitario en el nivel de Alcohol produce un *aumento del Beneficio esperado de 4.13 unidades**
Todas estas afirmaciones se hacen a constancia de todas las demás variables involucradas en el modelo (ceteris paribus), es decir, se valora el cambio marginal que tiene cada variable para valores fijos de todas las demás.
Así, ese aumento en los vinos de **** se da en aquellos que tienen la misma etiqueta, misma calificación del productor y mismos niveles de acidez y alcohol.